home *** CD-ROM | disk | FTP | other *** search
/ Transactor / Transactor_18_1987_Transactor_Publishing.d64 / n-body simulator (.txt) < prev    next >
Commodore BASIC  |  2023-02-26  |  8KB  |  319 lines
  1. 100 rem n-body simulator
  2. 110 rem version 6.09
  3. 120 rem by richard lucas
  4. 130 :
  5. 140 rem initialize
  6. 150 poke2038,peek(55):poke2039,peek(56)
  7. 160 poke 56,62: poke 55,0: clr
  8. 170 dim x(40),y(40),z(40),u(40),v(40),w(40),x1(40),y1(40),z1(40)
  9. 180 dim m(40),gm(40),e2(7)
  10. 190 dim x0(40),y0(40),z0(40),u0(40),v0(40),w0(40),a0(40),b0(40),c0(40)
  11. 200 dim ex(7),un(3),un$(3),tu$(3)
  12. 210 fori=1to3:readtu$(i),un$(i):next
  13. 220 data "days","10^7kilometers-tons-days","days","[193][213]-[197]arth mass-days","seconds"
  14. 230 data 1000km-kg-sec
  15. 240 gosub 2820
  16. 250 g=6.67e-11:sy=1:cf=40:c4=504:c7=7:c8=248:vi=53248:hi=vi+16:c5=255
  17. 260 hr=16*1024:o$="epncdls"+chr$(133)+chr$(136)+"xyzuvwmtr[196]":sp=0:dv=8
  18. 270 nb=0:dt=10:d2=dt*dt/2:d3=d2*dt/3:d4=d3*dt/4:cb=1
  19. 280 fori=0to7:ex(i)=2^(7-i):next
  20. 290 fori=0to7:e2(i)=2^i:next
  21. 300 un(1)=1000*86400^2/1e10^3
  22. 310 un(2)=5.9742e24*86400^2/1.4959789e11^3
  23. 320 un(3)=1/1e6^3
  24. 330 poke 53280,0:poke 53281,0:print chr$(14)+chr$(8)+chr$(151);
  25. 340 b$="":l$="":fori=1to38:b$=b$+chr$(32):l$=l$+chr$(192):nexti
  26. 350 b$=chr$(29)+b$+chr$(29)
  27. 360 :
  28. 370 rem main loop
  29. 380 gosub 2460
  30. 390 gosub 2620
  31. 400 print ">";
  32. 410 get a$:ifa$=""then 410
  33. 420 fori=1tolen(o$):ifa$=mid$(o$,i,1)then printa$:goto 450
  34. 430 next
  35. 440 goto 410
  36. 450 on i goto 490,510,1970,2900,3070,2110,2250,3190,3230
  37. 460 if i=18 then input"[206]ew name of body";n$(cb):goto 370
  38. 462 if i=19 then input"[196]rive #";dv:goto 370
  39. 470 input"new value";nv
  40. 480 on i-9 goto 2390,2400,2410,2420,2430,2440,2380,2450
  41. 490 poke55,peek(2038):poke56,peek(2039):clr:end
  42. 500 :
  43. 510 rem plot trajectories
  44. 520 if nb=0 then print "[206]o bodies in current system.":goto 370
  45. 530 if sp then 610:rem skip hires for sprites
  46. 540 rem set up hires screen
  47. 550 poke 56578,peek(56578)or3:rem switch to vic bank 1 (16k-32k)
  48. 560 poke 56576,(peek(56576)and252)or2
  49. 570 poke 53272,(peek(53272)and15)or128:rem char screen is in 9th k
  50. 580 poke 53265,peek(53265)or32:rem turn on hires screen
  51. 590 poke 820,0:poke 821,64:poke 822,0:poke 823,96:poke 251,0:sys 49152
  52. 600 poke 820,0:poke 821,96:poke 822,231:poke 823,99:poke 251,16:sys 49152
  53. 610 ifspthenfori=vitohi:pokei,.:next
  54. 620 ifspthen a1=0:fori=1to8:a1=a1 or(-(i<=nb)*2^(i-1)):next:poke vi+21,a1
  55. 630 ifspthenpoke 53281,0:printchr$(147)
  56. 640 t=0
  57. 650 :
  58. 660 rem move start parameters to working arrays
  59. 670 fori=1tonb
  60. 680 x(i)=x0(i)
  61. 690 y(i)=y0(i)
  62. 700 z(i)=z0(i)
  63. 710 u(i)=u0(i)
  64. 720 v(i)=v0(i)
  65. 730 w(i)=w0(i)
  66. 740 next
  67. 750 rem compute accel at time dt before start
  68. 760 fori=1tonb
  69. 770 a0(i)=.:b0(i)=.:c0(i)=.
  70. 780 next
  71. 790 fori=1tonb
  72. 800 ax=.:ay=.:az=.
  73. 810 forj=1tonb
  74. 820 ifi=jthen910
  75. 830 dx=x(j)-x(i)
  76. 840 dy=y(j)-y(i)
  77. 850 dz=z(j)-z(i)
  78. 860 r=sqr(dx*dx+dy*dy+dz*dz)
  79. 870 r3=r*r*r/gm(j)
  80. 880 ax=ax+dx/r3
  81. 890 ay=ay+dy/r3
  82. 900 az=az+dz/r3
  83. 910 next
  84. 920 x1(i)=x(i)-u(i)*dt+ax*d2
  85. 930 y1(i)=y(i)-v(i)*dt+ay*d2
  86. 940 z1(i)=z(i)-w(i)*dt+az*d2
  87. 950 next
  88. 960 fori=1tonb
  89. 970 ax=.:ay=.:az=.
  90. 980 forj=1tonb
  91. 990 ifi=jthen1080
  92. 1000 dx=x1(j)-x1(i)
  93. 1010 dy=y1(j)-y1(i)
  94. 1020 dz=z1(j)-z1(i)
  95. 1030 r=sqr(dx*dx+dy*dy+dz*dz)
  96. 1040 r3=r*r*r/gm(j)
  97. 1050 ax=ax+dx/r3
  98. 1060 ay=ay+dy/r3
  99. 1070 az=az+dz/r3
  100. 1080 next
  101. 1090 a0(i)=ax
  102. 1100 b0(i)=ay
  103. 1110 c0(i)=az
  104. 1120 next
  105. 1130 :
  106. 1140 rem calculate new system state
  107. 1150 fori=1tonb
  108. 1160 a1=.:b1=.:c1=.:a2=.:b2=.:c2=.
  109. 1170 forj=1tonb
  110. 1180 ifi=jthen1270
  111. 1190 dx=x(j)-x(i)
  112. 1200 dy=y(j)-y(i)
  113. 1210 dz=z(j)-z(i)
  114. 1220 r=sqr(dx*dx+dy*dy+dz*dz)
  115. 1230 r3=r*r*r/gm(j)
  116. 1240 a1=a1+dx/r3
  117. 1250 b1=b1+dy/r3
  118. 1260 c1=c1+dz/r3
  119. 1270 next
  120. 1280 j0=(a1-a0(i))/dt
  121. 1290 k0=(b1-b0(i))/dt
  122. 1300 l0=(c1-c0(i))/dt
  123. 1310 x2=x(i)+u(i)*dt+a1*d2+j0*d3
  124. 1320 y2=y(i)+v(i)*dt+b1*d2+k0*d3
  125. 1330 z2=z(i)+w(i)*dt+c1*d2+l0*d3
  126. 1340 forj=1tonb
  127. 1350 ifi=jthen1440
  128. 1360 dx=x(j)-x2
  129. 1370 dy=y(j)-y2
  130. 1380 dz=z(j)-z2
  131. 1390 r=sqr(dx*dx+dy*dy+dz*dz)
  132. 1400 r3=r*r*r/gm(j)
  133. 1410 a2=a2+dx/r3
  134. 1420 b2=b2+dy/r3
  135. 1430 c2=c2+dz/r3
  136. 1440 next
  137. 1450 j1=(a2-a1)/dt
  138. 1460 k1=(b2-b1)/dt
  139. 1470 l1=(c2-c1)/dt
  140. 1480 m1=(a2-2*a1+a0(i))/(dt*dt)
  141. 1490 n1=(b2-2*b1+b0(i))/(dt*dt)
  142. 1500 o1=(c2-2*c1+c0(i))/(dt*dt)
  143. 1510 x1(i)=x(i)+u(i)*dt+a1*d2+j1*d3+m1*d4
  144. 1520 y1(i)=y(i)+v(i)*dt+b1*d2+k1*d3+n1*d4
  145. 1530 z1(i)=z(i)+w(i)*dt+c1*d2+l1*d3+o1*d4
  146. 1540 u1(i)=u(i)+a1*dt+j1*d2+m1*d3
  147. 1550 v1(i)=v(i)+b1*dt+k1*d2+n1*d3
  148. 1560 w1(i)=w(i)+c1*dt+l1*d2+o1*d3
  149. 1570 a0(i)=a1
  150. 1580 b0(i)=b1
  151. 1590 c0(i)=c1
  152. 1600 next
  153. 1610 :
  154. 1620 fori=1tonb
  155. 1630 x(i)=x1(i)
  156. 1640 y(i)=y1(i)
  157. 1650 z(i)=z1(i)
  158. 1660 u(i)=u1(i)
  159. 1670 v(i)=v1(i)
  160. 1680 w(i)=w1(i)
  161. 1690 x=x(i):ifx(i)<.orx(i)>319then1720
  162. 1700 y=y(i):ify(i)<.ory(i)>199then1720
  163. 1710 gosub2710
  164. 1720 next
  165. 1730 t=t+dt
  166. 1740 if sp then printchr$(19);t
  167. 1750 geta$:ifa$=""then1150
  168. 1760 :
  169. 1770 rem restore character screen
  170. 1780 poke53265,peek(53265)and223
  171. 1790 poke 56578,peek(56578)or3
  172. 1800 poke 56576,(peek(56576)and252)or3
  173. 1810 poke 53272,(peek(53272)and15)or16
  174. 1820 poke vi+21,0:rem turn off sprites
  175. 1830 poke 53281,0
  176. 1840 if a$=chr$(133) then 1860
  177. 1850 goto 370
  178. 1860 print chr$(17);"[211]toring present system in memory.":gosub 3160
  179. 1870 for i=1 to nb
  180. 1880 x0(i)=x(i)
  181. 1890 y0(i)=y(i)
  182. 1900 z0(i)=z(i)
  183. 1910 u0(i)=u(i)
  184. 1920 v0(i)=v(i)
  185. 1930 w0(i)=w(i)
  186. 1940 next i
  187. 1950 goto 370
  188. 1960 :
  189. 1970 rem get new system from user
  190. 1980 input"[206]umber of bodies";nb
  191. 1990 ifnb<1ornb>50then370
  192. 2000 for i=1 to nb
  193. 2010 cb=i:gosub 2460
  194. 2020 print"[206]ame of body"i;:inputn$(i)
  195. 2030 iflen(n$(i))>25then2020
  196. 2040 print"[205]ass of body"i;
  197. 2050 input m(i):gm(i)=g*m(i)*un(sy)
  198. 2060 input"[201]nput location in x,y,z form";x0(i),y0(i),z0(i)
  199. 2070 input"[201]nput velocity in u,v,w form";u0(i),v0(i),w0(i)
  200. 2080 next i
  201. 2090 goto 370
  202. 2100 :
  203. 2110 rem load system description from disk
  204. 2120 print "[204]oad system data from disk."
  205. 2130 input"[212]ype name of data file";a$
  206. 2140 if len(a$)>13 then print"[212]oo long.":goto 2130
  207. 2150 open 15,dv,15:open 2,dv,2,"0:"+a$+".nb,s,r"
  208. 2160 gosub 3170
  209. 2170 if er<>0 then print er$(1);er$(2);er$(3);er$(4):gosub 3160:goto 2230
  210. 2180 input#2,sy,nb
  211. 2190 fori=1tonb
  212. 2200 input#2,n$(i),x0(i),y0(i),z0(i),u0(i),v0(i),w0(i),m(i)
  213. 2210 gm(i)=g*m(i)*un(sy)
  214. 2220 next
  215. 2230 close 2:close 15:cb=1
  216. 2240 goto 370
  217. 2250 rem save current system to disk
  218. 2260 print "[211]ave current system."
  219. 2270 input"[212]ype name of file";a$
  220. 2280 if len(a$)>13 then print"[206]ame too long.":goto 2250
  221. 2290 open 15,dv,15:c$=chr$(13)
  222. 2300 open 2,dv,2,"0:"+a$+".nb,s,w"
  223. 2310 gosub 3170
  224. 2320 if er<>0 then print er$(1);er$(2);er$(3);er$(4):gosub 3160:goto 2230
  225. 2330 print#2,sy;c$;nb
  226. 2340 fori=1tonb
  227. 2350 print#2,n$(i);c$;x0(i);c$;y0(i);c$;z0(i);c$;u0(i);c$;v0(i);c$;w0(i);c$;m(i)
  228. 2360 next
  229. 2370 close 2:close 15:goto 370
  230. 2380 m(cb)=nv:gm(cb)=g*m(cb)*un(sy):goto 370
  231. 2390 x0(cb)=nv:goto 370
  232. 2400 y0(cb)=nv:goto 370
  233. 2410 z0(cb)=nv:goto 370
  234. 2420 u0(cb)=nv:goto 370
  235. 2430 v0(cb)=nv:goto 370
  236. 2440 w0(cb)=nv:goto 370
  237. 2450 dt=nv:d2=dt*dt/2:d3=d2*dt/3:d4=d3*dt/4:goto 370
  238. 2460 rem display current system values
  239. 2470 printchr$(147)"            [206]-[194][207][196][217] [211][201][205][213][204][193][212][207][210]"
  240. 2480 printchr$(176);l$;chr$(174);
  241. 2490 l1$=chr$(221)
  242. 2500 printl1$" name: "n$(cb);tab(79);l1$;
  243. 2510 printl1$" body #"cb;tab(17)"mass:"m(cb);tab(39)l1$;
  244. 2520 printl1$" x:"x0(cb);tab(60)"u:"u0(cb);tab(79);l1$;
  245. 2530 printl1$" y:"y0(cb);tab(20)"v:"v0(cb);tab(39);l1$;
  246. 2540 printl1$" z:"z0(cb);tab(60)"w:"w0(cb);tab(79);l1$;
  247. 2550 printchr$(173);l$;chr$(189)
  248. 2560 print "number of bodies:"nb
  249. 2570 print "time interval:"dt;tu$(sy)
  250. 2580 print "unit system: "un$(sy)
  251. 2590 if sp then print "sprite mode";chr$(17):return
  252. 2600 print "hires point mode";chr$(17):return
  253. 2610 :
  254. 2620 rem display menu
  255. 2630 print "e[146]xit       p[146]lot       n[146]ew system"
  256. 2640 print "sc[146]ale      d[146]isplay    l[146]oad [196][146]rive#"
  257. 2650 print "s[146]ave       f1[146] prev body   f7[146] next body"
  258. 2660 print "x[146] position y[146] position z[146] position"
  259. 2670 print "u[146] velocity v[146] velocity w[146] velocity"
  260. 2680 print "m[146]ass       t[146]ime       r[146]ename"
  261. 2690 return
  262. 2700 rem plot point on hires screen
  263. 2710 if sp=1 then 2750
  264. 2720 ml=hr+(yandc8)*cf+(yandc7)+(xandc4)
  265. 2730 poke ml,peek(ml)orex(xandc7)
  266. 2740 return
  267. 2750 if i>8 then return
  268. 2760 x=x+24:y=y+50
  269. 2770 poke vi+(i-1)*2,xandc5
  270. 2780 poke vi+i*2-1,y
  271. 2790 ifx>c5thenpokehi,peek(hi)ore2(i-1)
  272. 2800 ifx<256thenpokehi,peek(hi)and(c5-e2(i-1))
  273. 2810 return
  274. 2820 rem ml code for high speed erase
  275. 2830 i=49152
  276. 2840 readmc:ifmc=256thenreturn
  277. 2850 pokei,mc:i=i+1:goto2840
  278. 2860 data173,52,3, 133,2, 173,53,3, 133,3
  279. 2870 data165,251, 160,0, 166,3
  280. 2880 data145,2, 236,55,3, 208,7, 166,2, 236,54,3, 240,9
  281. 2890 data230,2, 208,236, 230,3, 76,14,192, 96,  256
  282. 2900 print chr$(147);"[211]elect a system of units."
  283. 2910 print chr$(17);"1. 1 pixel = 10^7 kilometers"
  284. 2920 print "   1 mass = 1000 kilograms"
  285. 2930 print "   1 time = 1 day"
  286. 2940 print chr$(17);"2. 1 pixel = 1 [193][213] (earth radius)"
  287. 2950 print "   1 mass = 1 earth mass"
  288. 2960 print "   1 time = 1 day"
  289. 2970 print chr$(17);"3. 1 pixel = 1000 kilometers"
  290. 2980 print "   1 mass = 1 kilogram"
  291. 2990 print "   1 time = 1 second"
  292. 3000 print: input"[215]hich system";sy
  293. 3010 if sy<1orsy>3then370
  294. 3020 fori=1tonb
  295. 3030 gm(i)=g*m(i)*un(sy)
  296. 3040 next
  297. 3050 goto 370
  298. 3060 :
  299. 3070 rem switch plot systems
  300. 3080 if sp=1 then sp=0:goto 370
  301. 3090 sp=1
  302. 3100 fori=15872to15872+8*64:pokei,.:next:rem blank out sprite images
  303. 3110 fori=0to7:poke15872+i*64,224:poke15875+i*64,224:next:rem form dot shape
  304. 3120 fori=0to7:poke2040+i,248+i:next:rem set sprite pointers
  305. 3130 fori=0to7:pokevi+39+i,i+1:next:rem set sprite colors
  306. 3140 poke vi+29,0  :poke vi+23,0  :rem compress sprites
  307. 3150 goto 370
  308. 3160 forde=1to1500:next:return
  309. 3170 input#15,er$(1),er$(2),er$(3),er$(4)
  310. 3180 er=val(er$(1)):return
  311. 3190 cb=cb-1:if cb<1 then cb=nb
  312. 3200 printchr$(19)chr$(17)chr$(17);b$;b$;b$;b$;b$:printchr$(19)
  313. 3210 gosub 2480
  314. 3220 fori=1to6:printchr$(17);:nexti:printchr$(29);:goto410
  315. 3230 cb=cb+1:if cb>nb then cb=1
  316. 3240 goto 3200
  317. 3250 poke53272,peek(53272)and247
  318. 3260 poke53265,peek(53265)and223
  319.